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ABSTRACT 

Studies of the growth of cosmic perturbations are typically focused on galac- 
tic scales and above. In this paper we investigate the evolution of perturbations 
in baryons, photons, and dark matter for masses below 10 6 M Q (or wavenumbers 
above 100 Mpc -1 ). Fluctuations on these scales are of interest and importance 
because they grow to become the earliest collapsed objects and provide the first 
light sources in the so called dark ages. We investigate both the linear evolution 
and the second-order nonlinear effects arising from the coupling of large-scale 
velocity fields to small-scale perturbations in the baryon density and the electron 
ionization fraction. We find that this second order nonlinear coupling dominates 
the growth of perturbations with masses < 1O 3 M immediately after recombi- 
nation, enhancing the baryon fluctuation amplitudes by a factor of ~ 5, but the 
nonlinear effect does not persist at late times. 

1. Introduction 

Theories of structure formation are aimed at describing how gravity affects the distri- 
bution of matter in the universe. The currently accepted paradigm is that gravitational 
instabilities amplify fluctuations imprinted in the cosmological matter and radiation density 
in the early universe, and these fluctuations eventually evolve into the observed web of cos- 
mic structures. By assuming these fluctuations to be small, a linear perturbation theory can 
be formulated and solved. The linear theory for a perturbed Friedmann-Robertson- Walker 
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metric was developed by Lifshitz (1946). The theory and subsequent modifications have been 
applied to a large number of cosmological problems and are described in detail in several 
textbooks (e.g., Weinberg 1972; Peebles 1980; Peacock 1999). 

Most calculations and analyses based on the cosmological linear perturbation theory are 
focused on galactic scales or above at masses of M > 1O 9 M , or wavenumbers of k < 10 
Mpc -1 . This is likely due to the fact that quantities such as the mass fluctuation power 
spectrum can be determined observationally only on these large scales, and that the compu- 
tations become expensive at high wavenumbers due to the rapid oscillations in the perturbed 
fields. The theory, however, is equally valid for small scales, and the evolution of M < 10 6 M Q 
perturbations is key to the understanding of the formation of the first generation objects in 
the high redshift universe. A detailed study of perturbation modes up to k ~ 1000 Mpc^ 1 
has found interesting features in the linear baryon field (Yamamoto, Sugiyama, & Sato 1998). 
For example, the tight coupling between baryons and photons, which is assumed until the 
epoch of recombination in most calculations, breaks down well before recombination on these 
very small scales. This allows the baryon density field to grow substantially by recombination 
under the influence of the cold dark matter (CDM). This behavior is contrary to the common 
notion that the baryon fluctuations are wiped out by photon diffusion or Silk damping (Silk 
1968) on small scales. After recombination, on the other hand, perturbation modes with 
k > 300 Mpc -1 are below the Jeans mass and exhibit oscillatory behavior, in contrast to 
the lower-/c modes which grow via Jeans instability immediately after recombination. This 
retards the baryon growth until the mode becomes Jeans unstable again at a later time. 

Beyond the linear order, there have been recent discussions regarding nontrivial second- 
order effects that can enhance the linear growth of fluctuations on very small scales during 
recombination (Shaviv 1998; Liu et al. 2001). This nonlinear effect arises from the coupling 
of the large scale velocity fields to the perturbations in the baryon density 5b and the electron 
ionization fraction 5 Xe on small scales. An amplification of ~ 10 4 in 5b has been reported 
(Shaviv 1998), although a more in-depth analysis (Liu et al. 2001) has found that the 
dramatic increase in the baryon amplitude is a result of neglecting the electron-photon 
diffusion term, which when included, would reduce the nonlinear enhancement to a factor 
of ~ 10. More specifically, Shaviv (1998) showed that the time scale governing the growth 
of small scale perturbations in the baryon density is ~ Tr/to, where t r is the characteristic 
time scale for the second-order terms and r is the time scale for acoustic oscillations in 
the baryon fluid. He took the time scale td for momentum transfer between electrons and 
photons due to Thomson scattering to be much larger than r|,/r . Liu et al. (2001) pointed 
out that this assumption is invalid during recombination since Tp ~ t r/ t o and Td cannot 
be ignored. In addition, Shaviv (1998) used the Saha equation to evaluate 5 Xe which is 
needed to calculate tr. Liu et al. (2001) instead used a rate equation for 5 Xe because the 
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assumption of thermal equilibrium required for the validity of the Saha equation does not 
hold during recombination. They then assumed a wave form solution for the baryon density 
field 5b oc exp (out — ikx), and estimated the growth of 5b from the real part of u, which they 
obtained from the dispersion relation given by the baryon fluid equation. 

In this paper, we investigate this small-scale nonlinear phenomenon in its full extent 
by using the complete cosmological perturbation theory. Our work differs from Liu et al. 
(2001) in that instead of tracking only the real part of 5b from the dispersion relation, 
we calculate the evolution of the perturbations in all particle species (baryons, photons, 
CDM, and neutrinos) at high wavenumbers k > 1000 Mpc -1 by solving the complete set 
of Einstein and Boltzmann equations consistently. We use the linear Boltzmann code in 
the COSMIC package (Ma & Bertschinger 1995; Bertschinger 1995) and incorporate the 
aforementioned nonlinear effects by adding second-order contributions to the photon-electron 
momentum transfer terms due to Thomson scattering. For consistency, we include these 
terms in the evolution equations for photons as well as baryons. This approach allows us to 
take into account possible nonlinear coupling effects and to calculate the growth, decay, and 
oscillations in all particle species accurately. 

In § 2.1 we summarize the key portions of the linear cosmological perturbation theory. 
In § 2.2 we introduce the nonlinear coupling terms in the baryon and photon equations as 
well as the equations for the perturbed ionization fraction. Numerical results for the linear 
and the nonlinear evolution of four representative high-A; modes, k = 1000, 2500, 5000, and 
10000 Mpc^ 1 , are presented in § 3. The corresponding masses are ~ 1000, 64, 8, and 1 M & , 
respectively. § 4 includes a discussion of the physical origin of this phenomenon and the 
reason why a second order term could dominate the growth of perturbations at the era of 
recombination. 



2. Cosmological Perturbation Theory 

2.1. Linear Perturbations 

The discussion here is restricted to the scalar modes of metric perturbations. We will 
work in the conformal Newtonian gauge 

ds 2 = a 2 (r)[-(l + 2iP)dr 2 + (1 - 2(f))dx 2 ] , (1) 

where and ip are the two scalar potentials that characterize the perturbations, r is the 
conformal time, x are the spatial dimensions, and a is the scale factor. The scalar potential 
ip has a simple physical interpretation of being the gravitational potential in the Newtonian 



-4- 



limit, and <p ~ ^ in the absence of massive neutrinos. 

We solve the full set of linear evolution equations for small perturbations in the metric 
and the phase-space distributions of CDM, photons, baryons, and neutrinos. The full theory 
in this gauge is described in detail in Mukhanov et al. (1992) and Ma and Bertschinger 
(1995); here we write out only the lowest two velocity moments of the phase space distribution 
for CDM (subscript c) , baryons (b), and photons (7), which are the key quantities in this 
study. In /c-space, we have 

S c = —0 C + 30 , 

a 

8b = -6b + 30 , 

6\ = --6 b + c 2 s k 2 5 b + k 2 ^j + ^L a fi e a T (e^-e b ), (2) 
a 6pb 



4 

5 1 = --# 7 + 40, 



9 1 = k 2 (^5 y — <7 7 ) + k 2 ip + an e a T (6 b — 7 ) , 

where 5 is the density fluctuation, 6 = ikv is divergence of the fluid velocity, and <r 7 is the 
sheer stress of the photon field which is coupled to higher moments not written down here. 
On the right hand side, ot is the Thomson scattering cross section, c s is the baryon sound 
speed, n e is the mean electron number density, and p 7 and p b are the mean photon and 
baryon energy densities. Dots denote derivatives with respect to the conformal time, and 
the densities and wavenumbers are all comoving. Massive neutrinos have no direct coupling 
to the nonlinear terms considered here, so we do not include them. 

Additional equations are needed to calculate n e . We define the mean ionization fraction 

x e as 

_ n e n e TTip . . 

n H (l - y)pb 

where nn is the mean hydrogen number density, y is the primordial helium mass fraction, 
and m p is the proton mass. The value x e — 1 corresponds to complete hydrogen ionization 
and neutral helium, and x e can exceed unity when helium is ionized at an earlier time. We 
choose this convention for convenience because for the redshift of interest to this paper, 
hydrogen is the dominant source of free electrons. For times much before recombination, x e 
obeys the equilibrium Saha equation to a good approximation, 

_ J_ ( f^ekBTb \ 1 -13.6eV/k B T b 



-Y6.*5eV/k B T b U\ 

l-x e n H V 2tt/i 2 ' 
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Here T b is the baryon temperature, ks is the Boltzmann constant, m e is the electron mass, 
and h is the Planck's constant. During recombination the rapidly declining free electron den- 
sity leads to a breakdown of ionization equilibrium, and one must integrate the appropriate 
kinetic equations. We use the recombination rate equation (Peebles 1968; Spitzer 1978) 

ci IT 

— - = aC r [P(T b )(l - x e ) - n H a^ 2 \T b )x 2 e } , (5) 

GST 

where the collisional ionization rate from the ground state is 

P(T b ) = ( !! ^f?) 3/2 e^/K B T ba m {Tb) > (6) 

and the recombination rate to the excited states is 

(2) 647T e 4 / k B T b Y l/ \ /rr , , (T , n ,, 8 , fl^6eV\ 

Here e is the electron charge and c is the speed of light. The net recombination rate to the 
ground state is reduced by the fact that an atom in the n = 2 level may be ionized before it 
decays to the ground state. The Peebles' reduction factor, C r , is the ratio of the net decay 
rate to the sum of the decay and ionization rates from the n = 2 level, 

r = A Q + A 2s ^is / R s 

where 

^ 2 \T b )=f3(T b )e hc ^ T " x % A a = -^-, \ a = 1.216 x 10" 5 , (9) 

a 2 n ls X% 

A 2s ^i s is the rate at which net recombination occurs through two-photon decay from the 2s 
level and is equal to 8.227 s _1 , and n is is the number density of hydrogen atoms in the Is 
state. 



2.2. Second-Order Perturbations 

A complete list of second order nonlinear contributions in perturbation theory includes 
many terms, but the dominant contribution comes from the terms of order 5 b v b (e.g., Hu, 
Scott, and Silk 1994; Dodelson and Jubas 1995). The other second order terms are propor- 
tional to vl, which are suppressed because v b ~ S b /kr 5 b for the scales of interest here with 
kr ^> 1. It will be demonstrated later in the section that the perturbed ionization fraction 
5 Xe is of the same order of magnitude as 5 b for the scales of interest. We will therefore include 
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perturbations to both the ionization fraction and the baryon density as the terms giving rise 
to second order effects in the evolution equations. 

It is easier to determine such second-order terms in real space. The relevant first-order 
equations to be modified are the b and # 7 equations in (3), which in real space take the 
form 



v b = — v b - c s V5 b - V^> + — I an e cr T (t; 7 - v b ) , 
a 6p b 

v 7 = — V ~~ ~~ + an e a T (v b — v 7 ) . (10) 

Allowing for spatial fluctuations in the ionization fraction x e , 

x e = x e (l + 5 Xe ) , (11) 

the perturbations in the electron density then contain two terms arising from inhomogeneities 
in the baryon density and the ionization fraction, respectively: 

n e = n e (l + 5 ne ) = n e (l + 5 b + 5 Xe ) . (12) 

We obtain the second-order contributions to the Thomson scattering terms in equation (10) 
by replacing n e with n e , p b with p b = p b {l + 5 b ), and p 7 with p 7 = p 7 (l + 5 7 ): 

v b = --v b - c 2 s V8 b - VV> + ^z-an e a T (l + 5 Xe + <5 7 )(v 7 - v b ) , 
a op b 

t>' 7 = -V Q<5 7 - a^j - VV> + an e a T (l + S Xe + 5 b )(v b - u 7 ) . (13) 

We will drop the second-order term containing 5 1 in the v b equation because on small scales, 
5 7 <C 5 Xe ~ 5 b due to photon diffusion damping. Transforming back to k space, the nonlinear 
terms of 0(5v) become convolutions, and we have 

b = ~0b + c 2 s k 2 5 b + k 2 ij + ^an^T^ - b ) 
a 6p b 



3p& 



an e a T k J d 3 k' [v 7 (k') - v b (k')] 5 Xe (\k - k'\) 



# 7 = k 2 ( —Sj — o" 7 ) + k 2 tp + an e a T (0 b — 9 7 ) 



+an e a T k J d 3 k' [v b (k') - v 7 (k')] [S Xe (\k - k'\) + 5 b (\k - k'\)} , (14) 

We also need an evolution equation for the ionization fraction fluctuation 5 Xe . It can be 
obtained by perturbing the zeroth-order ionization rate equation (5) (Liu et al. 2001): 



dr 



-aC r [(35 Xe + n H a^x e (5 b + 25 X J] + - 5 x )j 1^ , (15) 
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where 

SC r A a ni s/ 9 (2) n H [(l - x e )5 b - x e 5 Xe 



C r [A a nis + A 2s - > i s n jf /(l - x e )] [K a rii s + (A 2s ^i s + /3( 2 ))n#(l - x e )] ' 



(16) 



We have assumed that 5x b = 0, i.e., there are no perturbations induced in the baryon 
temperature. We verify this claim by including the perturbed evolution equation for T b (Ma 
and Bertschinger 1995) in our calculations: 



d(6 Tb ) , 



T b _ 2d 
f b a 



+ ^^-an e a T [S Xe (T, - f b ) - T b 6 Tb ] . (17) 



Here /z is the mean molecular weight (including free electrons and all ions of Hydrogen and 
Helium) and T b = T b (l + JyJ.We find that the value for 5T b is always at least two orders 
of magnitude smaller than 5 Xe because the source term (T 7 — T b ) for ST b is negligible until 
well into the recombination era. Even after recombination the growth is not rapid enough 
to affect the evolution of 5 r . 

We note that the perturbation to the ionization fraction 8 Xe arises naturally from the 
first-order cosmological equations and is not an ad hoc term. This perturbation is ignored 
in standard linear calculations because its effects on the photon and baryon evolution enter 
only at second order, as shown by equations (13) and (14). We show in the next section that 
numerical results from the linear theory indeed give comparable amplitudes for 5 Xe and S b . 
We are thus justified in keeping both the 5 b v and 5 Xe v terms in our second order analysis. 



3. Numerical Results 

In this section we present results from numerical integration of the linear and second- 
order equations in § 2. The calculations are performed with the full Boltzmann code for the 
conformal Newtonian gauge in the COSMIC package (Ma & Bertschinger 1995; Bertschinger 
1995). For the second-order results, modifications are made to this code to include the 
nonlinear terms in equation (14) and the perturbed rate equation (15). The matter density 
parameter is taken to be fl m = 0.35 with Q c = 0.30 in CDM and Q b = 0.05 in baryons. The 
cosmological constant and Hubble parameter are Q\ = 0.65 and h = 0.75. Three species of 
massless neutrinos and a helium fraction y = 0.24 are assumed. The initial conditions are 
adiabatic, Gaussian fluctuations with a spectral index of n — 1. All results are normalized 
to COBE. 
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3.1. Linear Evolution 

Figure 1 illustrates the time evolution of the linear baryon density field 8b and CDM 
density field 8 C for wavenumbers k = 1000, 2500, 5000, and 10 4 Mpc -1 . The corresponding 
masses are ~ 1000, 64, 8, and 1 M , respectively. These scales are much smaller than the 
typically studied scales of k < 10 Mpc -1 or M > 1O 9 M . We note two interesting features 
that are absent in the more familiar behavior of 8b at lower k. First, even though photon 
diffusion damping wipes out 8b very rapidly at z ~ 10 6 to 10 5 , 8b at such high k is able 
to grow before recombination due to the break down of photon-electron tight coupling on 
small scales. This rapid regenerative growth is characterized by a (1 + z)~ 7 / 2 dependence 
and can be understood in terms of the balance achieved between the radiation drag force 
and the gravitational force (Yamamoto et al. 1998). As a result, 8b grows by at least 5 orders 
of magnitude before recombination ends at z ~ 1000. We also note a short period during 
recombination when 8b grows even more rapidly. This is because gravity is rapidly becoming 
the dominant force over the decreasing radiation drag force. 

The second feature special to these high k modes is that they undergo a second Jeans 
length crossing after recombination. It is well known that the baryonic Jeans wavenum- 
ber kj — Xj 1 ~ c s 1 [^nc?Gp) x l' 1 increases rapidly during recombination, reaching kj^ rec = 
900(fi m /i 2 ) 1 / 2 Mpc~ 1 ~ 370 Mpc -1 shortly after recombination for our model (Yamamoto et 
al. 1998). Perturbations with k > k Jrec , which were Jeans unstable before recombination, 
become stable after recombination and exhibit oscillations due to the baryon thermal pres- 
sure. The growth of 8b is therefore slowed down as can been seen in Figure 1. As we will show 
below, it is during this era when the small scale perturbations are relatively constant that 
the second order nonlinear term becomes important. When the Jeans wavenumber grows to 
the scale of a given k mode eventually, the growth of 8b picks up again. 

Figure 2 shows the linear power spectrum for the velocity difference (v b — i> 7 ) between 
the baryons and photons shortly after recombination. The Thomson scattering term is 
proportional to this quantity, which serves as the source term for the second-order effect to 
be discussed in the next section. 



3.2. Nonlinear Evolution 

We now solve the Einstein and Boltzmann equations taking into account the second- 
order terms in equation (14). Two new features must be handled: the evolution of the 
perturbed ionization fraction 8 Xe given by equation (15), and the /^-convolution in equa- 
tion (14). 
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To compute 5 Xe , we add equation (15) directly in the COSMIC Boltzmann code. We 
note that even though 5 Xe is a first-order quantity and is of comparable amplitude as 5b, 
equation (15) was not included in the original linear COSMIC code because 5 Xe contributes 
to the photon and baryon evolution equations only at second order. Figure 3 illustrates our 
numerical result for the time evolution of S Xe for the k = 5000 Mpc -1 mode. The dashed 
curve is for 5 Xe computed from the linear theory. It indeed has a similar amplitude as 5b in 
Figure 1. Including the second-order terms in the calculations enhance the amplitude of 5 Xe 
(solid curve), which peaks at a redshift of ~ 1000. 

To handle the convolution terms in equation (14), we use the property that the velocity 
difference (t> 7 — v^) in the integrand has significant contributions only from k' = 0.01 to 1 
Mpc" 1 , as shown in Figure 2. This is much smaller than k > 1000 Mpc" 1 investigated in 
this paper. This allows us to take out 5 Xe and 5b from the convolutions and approximate the 
velocity integral as 

J d 3 k' [v y (k') - v b (k')]6 Xe (\k - k'\) » a v 5 Xe (k) , (18) 

where 

a 2 v = 4n J dlnk'k' 3 [v 7 (k')-v b (k')} 2 . (19) 

Figure 4 shows a v as a function of redshift. Initially a v is small because the tight coupling 
between baryons and photons keeps v 1 ~ v b - As recombination proceeds and the tight 
coupling breaks down, the baryons fall into the CDM potential well, giving rise to the 
rapidly growing a v during z ~ 1500 to 1000. For comparison, the long-dashed curve shows 
the baryon velocity alone. The agreement between the two curves shows that the photon 
velocity field t> 7 is negligible after recombination. We approximate the time dependence of 
a v with the fitting formula 

log 10 a v = -3.57 + 9.56 x 10" 4 2 - 3.3 x 10~V + 1.126 x 10" V (20) 

shown as the short-dashed curve in Figure 4. This formula is used in the numerical integration 
for computational efficiency 

The results of our nonlinear calculations are illustrated in Figures 5-7. Figure 5 shows 
the time evolution of the linear vs. nonlinear baryon density field 5b for a single mode 
k = 5000 Mpc -1 ofMsi 1OM . The nonlinear 5 b is calculated from the full Boltzmann code 
with the second order convolution terms. The conformal time r is used for the horizontal 
axis to illustrate the post-recombination Jeans oscillations which are special to these high-A; 
modes (see also discussions in § 3.1). It can be seen that 5 b begins to oscillate immediately 
after recombination at r xs 200 Mpc due to the sharp decrease in the Jeans length. The 
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time period of these oscillations is roughly 2n/kc s , where the baryon sound speed is related 
to the baryon temperature by c 2 s = k B T b {l — d\nT b /d\na/3). For k = 5000 Mpc" 1 , the 
baryons have T b ~ 2500 K at time r ~ 200 Mpc, giving c s ~ 1.5 x 10~ 5 and At ~ 100 
Mpc, which agrees well with the periods of oscillations shown in Figure 5. Figure 6 shows 
the ratio of nonlinear and linear S b for four wavenumbers: k = 1000, 2500, 5000, and 10000 
Mpc -1 . Figure 7 shows the same thing for the photon density field <5 7 , which has a similar 
behavior as the baryons. 

We find that the inclusion of the second-order terms has two main effects on the baryons: 
an increase in the amplitude of the Jeans oscillations, and a slight shift in the sound speed 
and hence the oscillation period. Both effects are clearly seen in the figures. The nonlin- 
ear enhancement reaches a maximum of a factor of several right after recombination. It 
then decreases as the mean ionization fluctuation x e decreases. The slight shift in the os- 
cillation period is caused by the inclusion of perturbations to the ionization fraction. We 
have effectively introduced another term in the evolution equation for 6 b , which together 
with the linear terms determine the period of the oscillations. The nonlinear coupling thus 
changes the time period of these oscillations to ~ 2-n-/(kc s + nonlinear effects). This slight 
shift in the effective sound speed is the reason for the oscillations in Figures 6 and 7. It is 
also worthwhile to note that the maximum amplitude achieved by the ratio of the baryonic 
perturbations in the nonlinear theory compared to the linear calculation does not follow a 
simple relationship as a function of k. The amplitude of the ratio, for example, is larger 
for k = 5000 Mpc -1 than for k = 10000 Mpc -1 . This is because the effect of the nonlinear 
term becomes less important as the frequency of the Jeans oscillations increases. For the 
particular cosmological model under consideration this occurs at a scale between k = 5000 
and 10000 Mpc -1 . 



4. Discussion 

Figures 6 and 7 show that the second order terms arising from the coupling of the large 
scale baryon velocity fields to the perturbation in the electron ionization fraction on small 
scales dominate the evolution equation for the baryon density perturbations 5 b at z ~ 1000. 
This second order effect dies off at lower z, however, and it does not push the small scale 
baryon fluctuations into the nonlinear regime at a significantly earlier time as was suggested 
by Shaviv (1998). 

It is nevertheless interesting that for modes with k > 1000 Mpc" 1 , a second-order term 
can significantly alter the growth predicted by the linear theory right after recombination. 
How does any second-order term at z ~ 1000 lead to an amplification of a factor of ~ 5 in 
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the linear baryon fluctuations as we have found? To understand this, we plot in Figure 8 the 
time evolution of all the terms on the right-hand-side of the 9b equation from (14): —a9b/a 
(short-dashed), c 2 k 2 Sb (dotted), (4p 1 /3p b )an e aT(9 1 — 9b) (long dashed), and the second-order 
term ~ VbS Xe (dot-dashed). Since the potential term k 2 i(j is nearly a constant of time, it is 
convenient to normalize other terms to k 2 ip. To contrast the scale dependence of the nonlinear 
effect, we show two modes: k = 250 and 5000 Mpc -1 . As one can see, the nonlinear effect is 
important only for the higher-/c mode. The nonlinear term is unimportant for the k = 250 
Mpc -1 mode because the pressure term k 2 c 2 Sb is small on these scales and can not prevent the 
perturbations from growing under the influence of gravity. For the k = 5000 Mpc -1 mode, 
however, the pressure term k 2 c 2 5b keeps the perturbations in the terminal velocity stage for 
a longer period of time (see Yamamoto et al. 1998). This leads to a near cancellation of all 
linear terms governing the growth of 9b- As a result, the second-order term gains importance 
and through the nonlinear coupling becomes the dominant term shortly after recombination. 
The effect of the coupling on the perturbed ionization fraction is shown in Figure 3, where 
the nonlinear coupling can give rise to an order of magnitude change in the value of 5 Xe . At 
late times the second-order term becomes negligible as the free electron density n e becomes 
very small. We note that the Hubble drag term —a9b/a is unimportant throughout because 
the scales of the perturbations under consideration are so small that they are not influenced 
by the large scale expansion of the universe. 

Our approach and results differ from the previous work on this subject (Liu et al. 2001) 
in a number of ways. We have solved the full set of Einstein and Boltzmann equations for the 
metric perturbations and the perturbation field in CDM, baryons, photons, and neutrinos. 
This approach has allowed us to calculate accurately the phase space distributions of both 
baryons and photons. For the baryons, we find that the maximum amplitude reached by 
the nonlinear to linear ratio, ^(Nonlinear) /^(Linear), is similar to that of Liu et al. (2001) 
for k < 5000 Mpc -1 , but our results show an oscillating ratio because the addition of the 
nonlinear coupling causes a change in the effective speed of sound. This effect was unnoticed 
in the work of Liu et al. (2001). In addition, their analysis finds a monotonic increase in 
the amplitude of the ratio as a function of k. We find instead that the maximum amplitude 
reached by the ratio for k = 10000 Mpc -1 is less than that for k = 5000 Mpc -1 . This is 
again because their analysis was unable to track the Jeans oscillations after recombination 
accurately. 

Our calculation has also included the nonlinear coupling terms in the evolution of the 
photon density perturbations. As equations (13) and (14) show, two comparable nonlinear 
terms contribute to the photon velocity field: 5 Xf ,{vb — v-y) and Sb(vb — t> 7 ). We have shown 
that the 5 Xe (vb — 1> 7 ) term results in an initial enhancement of the amplitude for the baryons. 
Momentum conservation in Thomson scattering implies that this nonlinear term would tend 
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to suppress the amplitude obtained from the linear calculation for the photons. It turns 
out that the 5b(vb — f 7 ) term is generally larger and it acts in an opposite sense as the 
$x e ( v b — Vy) term. The end result is an enhancement in <5 7 as well, as shown in Figure 7. It 
may be an interesting future study to compute directly the angular power spectrum for the 
temperature fluctuations in the cosmic microwave background, but we expect the angular 
scales associated with the k > 1000 Mpc -1 modes to be at multipoles of I > 10 6 , below the 
scales probed by the current generation of experiments. 

We thank Ue-Li Pen, Naoshi Sugiyama, and Thanu Padmanabhan for useful discussions. 
C.-P. M. acknowledges support of an Alfred P. Sloan Foundation Fellowship, a Cottrell 
Scholars Award from the Research Corporation, a Penn Research Foundation Award, and 
NSF grant AST 9973461. 
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Fig. 1. — Time evolution of the linear baryon (lower curve) and cold dark matter (upper 
curve) density field on very small scales: k = 1000, 2500, 5000, and 10000 Mpc -1 . The results 
are computed in the conformal Newtonian gauge using the COSMIC Boltzmann code. The 
cosmological parameters are Q c = 0.30, fib = 0.05, = 0.65, h = 0.75, and the fluctuation 
amplitudes are normalized to COBE. Photon diffusion damping results in the sharp drop in 
the baryon amplitude at z ~ 10 6 to 10 5 . This is followed by a power law growth when the 
photon-baryon tight coupling becomes ineffective. The small oscillations after recombination 
at z ~ 1000 in the high-fc modes are a result of the second Jeans length crossing (see text). 
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Fig. 2. — Linear power spectrum for the velocity difference between baryons and photons at 
redshift z = 800. The dominant contributions come from the modes in the range of k ~ 0.01 
to 1 Mpc -1 , which is much below the modes of interest in Fig. 1. This feature allows us to 
simplify the second-order velocity convolution terms in eq. (14). 
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Fig. 3. — Perturbation to the ionization fraction, 5 Xe , as a function of redshift, computed 
without (dashed curve) and with (solid curve) the second-order terms in eq. (14) for the 
velocity field 9. The peak occurs z ~ 1000, which is the same for all k modes investigated. 
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Fig. 4. — Evolution of a v , the rms of the baryon velocity (long dashed curve) and the 
rms of the baryon and photon velocity difference (vb — f 7 ) (solid curve). The velocity dif- 
ference is small at z > 1500 due to the tight coupling but it increases rapidly during and 
after recombination when baryons and photons are no longer tightly coupled. The baryons 
subsequently fall into the potential well of the CDM and results in an increase in Vb and 
(vb — Vy), whereas the photon perturbations simply diffuse away. The short-dashed curve 
shows the fitting formula for fob — u 7 ) used in our calculations. 
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Fig. 5. — The linear (solid curve) vs. nonlinear (dashed curve) baryon density field 8^ as 
a function of the conformal time r for the k = 5000 Mpc -1 mode. The Jeans oscillations 
discussed in the text are clearly seen. The effect of the second-order terms is to increase 
the oscillation amplitude as well as shift slightly the position of the peaks. The oscillation 
period is ~ 90 Mpc immediately after recombination. The conformal time is r w 200 Mpc 
at recombination and r ~ 10 4 Mpc today. 
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Fig. 6. — Ratio of the nonlinear and linear baryon density field 8b for four k modes: 1000, 
2500, 5000, and 10000 Mpc" 1 . The oscillations in the ratio are more rapid for higher k 
modes due to the faster Jeans oscillations. The maximum amplitude reached by a particular 
k mode is not a monotonic function of k but peaks between k = 5000 and 10000 Mpc -1 . 
The oscillations die out and the effect of the nonlinear term decreases at low redshifts when 
other terms begin to dominate. 
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Fig. 7. — Same as Fig. 6 but for the photon density field 5 1 . The behavior is qualitatively 
similar to the baryons in Fig. 6. 
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Fig. 8. — Amplitude of the individual terms on the right-hand side of the baryon velocity 
equation (14) relative to k 2 ip: —a6b/a (short-dashed), c 2 k 2 Sb (dotted), (4p 7 /3pb)ara e o"r(^ 7 — 
6b) (long dashed), and the second-order term (dot-dashed). For the lower k mode, the 
second-order term is never important compared with the first-order terms. For k = 5000 
Mpc -1 , however, the near cancellation among the first-order terms at z ~ 1000 allows the 
second order term to dominate, leading to the nonlinear enhancement shown in Figs. 6 and 
7. This figure illustrates how the second order term becomes important for high k modes. 



